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We consider a class of nonlinear stochastic differential equations, giving the power-law behavior of 
the power spectral density in any desirably wide range of frequency. Such equations were obtained 
starting from the point process models of 1/f 13 noise. In this article the power-law behavior of 
spectrum is derived directly from the stochastic differential equations, without using the point 
process models. The analysis reveals that the power spectrum may be represented as a sum of the 
Lorentzian spectra. Such a derivation provides additional justification of equations, expands the 
class of equations generating 1/ f 13 noise, and provides further insights into the origin of l//' 3 noise. 
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I. INTRODUCTION. 



Power-law distributions of spectra of signals, including 1/f noise (also known as 1// fluctuations, nicker noise 
and pink noise), as well as scaling behavior in general, are ubiquitous in physics and in many other fields, including 
natural phenomena, human activities, traffics in computer networks and financial markets. This subject has been a hot 
research topic for many decades (see, e.g., a bibliographic list of papers by Li [l[, and a short review in Scholarpedia 

• : !)• 

Despite the numerous models and theories proposed since its discovery more than 80 years ago .3, 4], the intrinsic 
origin of 1// noise still remains an open question. There is no conventional picture of the phenomenon and the 
i " i ' mechanism leading to 1/f fluctuations are not often clear. Most of the models and theories have restricted validity 
because of the assumptions specific to the problem under consideration. A short categorization of the theories and 
t— I ■ models of 1/f noise is presented in the introduction of the paper [j| . 

Until recently, probably the most general and common models, theories and explanations of 1// noise have been 
based on some formal mathematical description such as fractional Brownian motion, the half-integral of the white 
noise, or some algorithms for generation of signals with scaled properties [fj|-fl4j and the popular modeling of 1/f 
noise as the superposition of independent elementary processes with the Lorentzian spectra and a proper distribution 
of relaxation times, e.g., a 1/T re i a x distribution [15l - l21| . The weakness of the latter approach is that the simulation 
of 1//^ noise with the desirable slope (3 requires finding the special distributions of parameters of the system under 
consideration; at least a wide range of relaxation time constants should be assumed in order to allow correlation with 
experiments |2a - l28j . 

Nonlinear stochastic differential equation with linear noise and non-linear drift, was considered in Ref. [9J. It was 
found that if the damping is decreasing with increase of the absolute value of the stochastic variable, then the solution 
of such a nonlinear SDE has long correlation time. Recently nonlinear SDEs generating signals with 1/f noise were 
obtained in Refs. [29|,[3(| (see also recent papers |l,|3l|), starting from the point process model of 1/f noise [27l I32rl39| . 

The purpose of this article is to derive the behavior of the power spectral density directly from the SDE, without 
using the point process model. Such a derivation offers additional justification of the proposed SDE and provides 
further insights into the origin of 1// noise. 
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II. PROPOSED STOCHASTIC DIFFERENTIAL EQUATIONS 



Starting from the point process model, proposed and analyzed in Refs. [271 |32H38| . the nonlinear stochastic differ- 
ential equations are derived 0, H^, [3(| ■ The general expression for the SDE is 



dx = 



a 2 (ri-~ S )x 2r '- 1 dt + ax r >dW . (1) 



'Electronic address: julius.ruscckas@tfai.vu.lt 



2 



Here x is the signal, rj is the exponent of the multiplicative noise, v defines the behavior of stationary probability 
distribution, and W is a standard Wiener process. 

SDE ([1]) has the simplest form of the multiplicative noise term, ax v dW. Multiplicative equations with the drift coef- 
ficient proportional to the Stratonovich drift correction for transformation from the Stratonovich to the Ito stochastic 
equation [40( generate signals with the power-law distributions |5| . Eq. (JlJ is of such type and has probability distribu- 
tion of the power-law form P(x) ~ x~ u . Because of the divergence of the power-law distribution and the requirement 
of the stationarity of the process, the SDE (fTJ) should be analyzed together with the appropriate restrictions of the 
diffusion in some finite interval. For simplicity, in this article we will adopt reflective boundary conditions at x = x m - m 
and x = x max . However, other forms of restrictions are possible. For example, exponential restriction of the diffusion 
can be obtained by introducing additional terms in Eq. (fTJ), 

dx = a 2 (V^ + ?(— I"™? ( —Y)* b '- 1 dt + *x'dW. (2) 



2 2 \ x J 2 
Here m is some parameter. 

Equation (fTJ) with the reflective boundary condition at x m i n and x max can be rewritten in a form that does not 
contain parameters a and x m i n . Introducing the scaled stochastic variable x — > x/x m - m and scaled time t — > a 2 x^^ t 
one transforms Eq. (JTJ to 

dx= (rj- ~) x^dt + x^dW . (3) 
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The scaled equation (J3)) has a boundary at x = 1 and at 

£ = — . (4) 

^min 

Further we will consider Eq. ([3]) only. In order to obtain noise we require that the region of diffusion of the 

stochastic variable x should be large. Therefore, we assume that £ 3> 1. 



III. POWER SPECTRAL DENSITY FROM THE FOKKER-PLANCK EQUATION 

According to Wiener-Khintchine relations, the power spectral density is 

/oo poo 
C(t)e lult dt = 4 C(t)cos(ut)dt, (5) 
-oo JO 

where uj = 2nf and C(t) is the autocorrelation function. For the stationary process the autocorrelation function can 
be expressed as an average over realizations of the stochastic process, 

C(t) = (x(t')x(t' + *)) . (6) 

This average can be written as 

C{t)= ( dx { dx'xx'P (x)P x (x',t\x,0), (7) 



where Pq(x) is the steady-state probability distribution function and P x {x' , t\x, 0) is the transition probability (the 
conditional probability that at time t the signal has value x' with the condition that at time t = the signal had the 
value x). The transition probability can be obtained from the solution of the Fokker-Planck equation with the initial 
condition P x (x' , 0\x, 0) = 6(x' — x). 

Therefore, for the calculation of the power spectral density of the signal x we will use the Fokker-Planck equation 
instead of stochastic differential equation ([3]) . The Fokker-Planck equation corresponding to the Ito solution of Eq. ([3]) 

is mm 

^P=-(v-^)^- l P+\?-^ r >P- (8) 
dt V ' 2 / dx 2dx 2 w 



The steady-state solution of Eq. ([8]) has the form 
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The boundary conditions for Eq. (|SJ) can be expressed using the probability current [42[ 

S(x,t) = („- ~) x 2 "^P- \^ 2r 'P- (10) 

At the reflective boundaries x m i n = 1 and a; max = £ the probability current S(x, t) should vanish, and, therefore, the 
boundary conditions for Eq. ([8]) are 

S(l,t) = 0, S(£,t) = 0. (11) 



A. Eigenfunction expansion 

We solve Eq. © using the method of eigenf unctions. An ansatz of the form 

P(x,t) =P x (x)e- Xf (12) 

leads to the equation 



('-5) 



where P\(x) are the eigenfunctions and A > are the corresponding eigenvalues. The eigenfunctions P\(x) obey the 
orthonormality relation [42j 



S 

e*MPx(x)Px,(x)dx = 5 x ,v, (14) 



'1 

where Q(x) is the potential, associated with Eq. ((HJ, 

<!>(x) = -ln.Po(x). (15) 

It should be noted that the restriction of diffusion of the variable x by x m - m and x max ensures that the eigenvalue 
spectrum is discrete. Expansion of the transition probability density in a series of the eigenfunctions has the form [42| 

P x (x',t\x,0) =Y,P^')^ ix) P^)^ Xt ■ (16) 

A 

Substituting Eq. (fTB]) into Eq. (J7J) we get the autocorrelation function 

C(t)=]T e - At A A 2 . (17) 

A 

Here 

X x = J xP x {x)dx (18) 

is the first moment of the stochastic variable x evaluated jwith the A-th eigenfunction P\(x). Such an expression for 
the autocorrelation function has been obtained in Ref. [43j. Using Eqs. (JSJl and (fT7| we obtain the power spectral 
density 

This expression for the power spectral density resembles the models of 1/f noise using the sum of the Lorentzian 
spectra (l5l - [T8l l27l . l28l . l44l |45| . Here we see that the Lorentzians can arise from the single nonlinear stochastic 
differential equation. Similar expression for the spectrum has been obtained in Ref. [l4| where reversible Markov 
chains on finite state spaces were considered (Eq. (34) in Ref. [lH with — 7fe, m playing the role of A ). 

A pure l//' 3 power spectrum is physically impossible because the total power would be infinity. It should be 
noted that the spectrum of signal x, obeying SDE ([3]), has 1//^ behavior only in some intermediate region of 
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frequencies, f m - ln <C / <C /max, whereas for small frequencies / <C /mm the spectrum is bounded. The behavior of 
spectrum at frequencies / m j n <C / <C /max is connected with the behavior of the autocorrelation function at times 
V/max <C t <C l//min- Often l//' 3 noise is described by a long-memory process, characterized by S(f) ~ 1//^ 
as / — > 0. An Abelian-Tauberian theorem relating regularly varying tails shows that this long-range dependence 
property is equivalent to similar behavior of autocorrelation function C(t) as t — > oo [46]. However, this behavior of 
the autocorrelation function is not necessary for obtaining required form of the power spectrum in a finite interval of 
the frequencies which does not include zero (47l - l49| . 

From Eq. (|19p it follows that if the terms with small A dominate the sum, then one obtains l// 2 behavior of the 
spectrum for large frequencies /. If the terms with 1//^ (with < 2) are present in Eq. (fT9|) . then at sufficiently large 
frequencies those terms will dominate over the terms with 1/f 2 behavior. Since the terms with small A lead to l// 2 
behavior of the spectrum, we can expect to obtain 1/ f& spectum in a frequency region where the main contribution 
to the sum in Eqs. (fTT|) and (flUf is from the large values of A. Thus we need to determine the behavior of the 
eigenfunctions P\{x) for large A. The conditions when eigenvalue A can be considered as large will be investigated 
below. 



B. Eigenfunctions of the Fokker-Planck equation 



For r\ ^ 1, it is convenient to solve Eq. (1131) by writing the eigenfunctions P\(x) in the form 

P x {x)=x- v u x (x x -v). (20) 
The functions u\[z) with z = x 1 ^ n obey the equation 

-^ux(z) - (2a - l)~ Ux (z) = -p 2 u x (z) , (21) 

where the coefficients a and p are 



v - 1 V2A , x 

° = 1 + 2(T^' " = iiFTJ- < 22 > 

The area of diffusion of the variable z = x l ~ r> is restricted by the minimum and maximum values z min and z max , 

Zmin_ \i, v <i, Zmax ~\e-\ »?<i. [ ' 

The probability current S x (x), Eq. (|10j) . rewritten in terms of functions u x , is 

S x {z) = \{ V -l)z^^-u x {z). (24) 

Therefore, the boundary conditions for Eq. (f2"Tj) . according to Eci. pT|) are tt A (l) = and "^(C 1 '') = 0- Here u' x (z) is 
the derivative of the function u x (z). The orthonormality relation (|14[) yields the orthonormality relation for functions 

ma(^), 

z t =*u x (z)u X r{z)dz = S Xi y . (25) 



The expression (|18|) for the first moment X x of the stochastic variable x evaluated with the A-th eigenfunction becomes 

1 f^ v i+ 

X X = - / z^Tu x (z)dz. (26) 



C. Solution of the equation for eigenfunctions 

The solutions of equation (|2TT) are [13] 

ua(z) = z a [dJ a (pz) + c 2 Y a (pz)} , (27) 
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where J a (z) and Y a (z) are the Bessel functions of the first and second kind, respectively. The coefficients c\ and 
C2 needs to be determined from the boundary and normalization conditions for function u\(z). The asymptotic 
expression for the function u\(z) is 

u\(z) w c\z a ~^p~^ cos(pz + a) , pz > 1 . (28) 

Here a is a constant to be determined from the boundary conditions and c\ is the constant to be determined from 
the normalization (|23|) . 

The behavior of the power spectral density in Eq. (fT9|) as 1//^ can be only due to terms with large A. Therefore, we 
will consider the values of A for which at least the product pz m ax is large, pz max 3> 1. The first moment of the variable 
x in the expression for the autocorrelation function (fTTj) is expressed via integral (f2l>|) . If the condition pz 3> 1 is 
satisfied for all z then the function u\(z) has frequent oscillations in all the region of the integration, and the integral 
is almost zero. Consequently, the biggest contribution to the sum in Eq. (|17[) makes the terms corresponding to those 
values of A, for which the condition pz 3> 1 is not satisfied for all values of z between z m i n and z max - Therefore, we 
will restrict the values of A by the condition pz mm <C 1. Thus we will consider eigenvalues A satisfying the conditions 

1 Amax < P « 1/Zmin ■ (29) 

Explicitly, we have conditions l<p< £»7 -1 if 17 > 1 and <C /? -C 1 if r\ < 1. 

The derivative of the function ua(z), Eq. (|27|) . is 

u a(*0 = [ciJa-i(pz) + c 2 F a _i(/jz)] . (30) 

Since we consider the case pz min <C 1, then, using Eq. ([3U)l . instead of the boundary condition w A (z m ; n ) = we can 
approximately take the condition 

lim (aJa-xiy) +c 2 y o -i(y)) =0. 

If a > 1 then we get C2 = 0; if a < 1 then C2 = — c\ tan(7ra). Using those values of the coefficient C2 we obtain the 
solutions of Eq. (|2"T|) 

(\ ~ J c\z a J- a (pz) , a<l, 
^c A z a J Q (pz) , a > 1 . 

From approximate solution (|31j) . using asymptotic expression for the Bessel functions, we can determine the parameter 
a in the asymptotic expression (|2"8")) . We obtain that the parameter a depnends on a and does not depend on p. 

D. Normalization 

Taking A = A' from Eq. (|25|) we get the normalization condition. Using Eq. (l3lT) we have 



zJ± a (pz)dz = — - — ^— - / yJl a (y)dy « 1 



(„-!)(!-,,) ^ ' („-l) |l-,|p 



Taking into account the condition pz m - m <C 1 and replacing the lower limit of integration by 0, we obtain that the 
integral is approximately equal to 

p2max 

yJ± a (y) d y 







7T 

Here we assumed that pz max 3> 1. Therefore, the normalization constant c\ is 




(32) 
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IV. CALCULATION OF THE POWER SPECTRAL DENSITY 
A. Estimation of the first moment X\ of the stochastic variable x 

The expression (fT7| for the autocorrelation function contains the first moment X\ of the variable x, expressed as 
an integral of the function u\, Eq. ([251) . Using Eq. (j3"Tl we get 

X\f=i- / z p ~ 1 J± a (pz)dz = - — — — =■ / y l3 ~ 1 J± a (y)dy . 



Here 

R = 1 -I- 

2(t? -1) 



/3 = 1 + 7T^4t (33) 



and "+" sign is for a > 1, while "— " is for a < 1. 

If ±q + /? > 0, taking into account that pz m - m <C 1, we can replace the lower limit of the integration by 0, 
X\ R3 -pr J PZmax y /3 ~ 1 J± a {y)dy . We get that the integral in the expression for X\ approximately does not 
depend on the lower limit of integration pz m in- We can integrate the integral by parts and use the properties of the 
Bessel functions to obtain 

Xx « {1 (t y fl - 1 J±( a -i ) (y)QZ ± 03 + « - 2) /" 2 J ±( a-i)(y)rfy) • (34) 

Using expression <f3T|) for the function -ma (2), the boundary conditions u A (z m in) = and it A (z m ax) = leads to 

J±(a-l){pZmm) = , J±(a-1) (/O^max) = . (35) 

Therefore, the first term in the expression (|34l) for X> is zero. If f3 < |, taking into account that pz max 3> 1, we can 
extend the upper limit of integration to +00. We get that the integral for X\ approximately does not depend on the 
upper limit of integration pz max . 

Therefore, the first moment X\ of the variable x is proportional to the expression 

|l-7/| pP 

Now we are ready to estimate the power spectral density. 



B. Power spectral density 

Since pz max 1 , from the boundary condition u' x {z mayi ) = using the asymptotic expression (|28p for the function 
u\{z) we obtain the condition sin(pz max + a) = and pz n 



A„«^-^(™-a) 2 . (37) 
zz. 



max 



Equation (j37j) shows that the density of eigvalues D(X) is proportional to 1/vA. Since the parameter o does not 
depend on A, it follows that the density of eigenvalues and, consequently, the autocorrelation function do not depend 
on th parameter a. 

In order to estimate the sum in expression (j!7[) for the autocorrelation function, we replace summation by the 
integration, 

C(t) w / e- xt XlD{X)d\ (38) 



Such a replacement is valid when pz max ^> 1. Using the approximate expressions (|32|) and (|36|) we get the expression 
for the autocorrelation function 

— 2 

C(t) ~ / 2mi " \~P e~ xt d\ = t?- 1 [r(l - /3, z-^i) - T(l - /3, z-f n i)] . (39) 
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Here T(a, z) = J z °° t a 1 e t dt is the incomplete Gamma function. When z^ in <f< we have the following lowest 
powers in the expansion of the approximate expression (I39p for the autocorrelation function in the power series of t: 



C(t) 



0-1 13-2 



/3 > 2 



+ (7 -l)t + t Intent), = 2 
^+^-^(1-/3), l</3<2 (40) 

l* r(i-0), /3<i 



Here 7 sa 0.577216 is the Euler's constant. Similar first terms in the expansion of the autocorrelation function in the 
power series of time t has been obtained in Rcf . Q . 

Similarly, when /9z max ^> 1, replacing in Eq. (|19[) the summation by the integration we obtain the power spectral 
density 



5(/)«4y ^—pXiD(\)d\. (41) 

Equation, similar to Eg. (|41j) has been obtained in Ref. @ by considering a relaxing linear system driven by white 
noise (Eq. (27) in Ref. Q). Similar equation also has been obtained in Ref. [l4| where reversible Markov chains on 
finite state spaces were considered. In both Refs. ||, [3] the power spectral density is expressed as a sum or an 
integral over the eigenvalues of a matrix describing transitions in the system. 
Using the approximate expressions (1321) and (|36[) we get the equation 

-2 

sif ^lZ^^ dX - (42) 

When z~^ x «w« z~ 2 n then the leading term in the expansion of the approximate expression (|4"2")l for the power 
spectral density in the power series of u is 

( , ,-13 a ^ 1 

(43) 

The second term in the expansion is proportional to oj~ 2 . In the case of (3 < 2, the term with oj~@ becomes larger than 
the term with lu~ 2 when z~ax ^ w - Therefore, we obtain 1//^ spectrum in the frequency interval 1 « w « ^C 7 ? -1 ) 
if 77 > 1 and the frequency interval l/£ 2 ( 1- '') <c u) <C 1 if r\ < 1. It should be noted that time t and frequency ui in 
our analysis are dimensionless. 

Eq. (|4ip shows that the shape of the power spectrum depends on the behavior of the eigenfunctions and the 
eigenvalues in terms of the function X 2 D(X). This function X\D{\) should be proportional to A - * 3 in order to obtain 
1//^ behavior. Similar condition has been obtained in Ref. 1141. E q. (13) for discrete time process and the unnumbered 
equation after Eq. (34) for a continuous time process in Ref. [14j are analogous to the condition XfD(X) ~ A _/3 since 
the density of eigenvalues in Ref. [13] is proportional to l/|7'(x)|. Equations ([39]) and (|42|) for the power spectral 
density and autocorrelation function resembles those obtained from the sum of Lorentzian signals with appropriate 
weights in Ref. [27j]. 




V. NUMERICAL EXAMPLES 



If v = 3 we get that (3 — 1 and stochastic differential equation ([T]) should give signal exhibiting 1// noise. We 
will solve numerically two cases: 77 = | > 1 and rj = — i < 1. For the numerical solution we use Euler-Marujama 
approximation, transforming differential equations to difference equations. Equation (1441) with r\ = 5/2 was solved 
using variable step of integration, solution Eq. (|45p with 77 = —1/2 was performed using a fixed step of integration. 

When T) = I and v = 3 then equation ([1]) is dx = a 2 x A dt + axidW . Using exponential restriction of the diffusion 
region we have the equation 

dx = <7 2 (l + l^L-^—^xUt + axUw. (44) 

\ 2 X 2 Xmax / 
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FIG. 1: (Color online) Probability distribution function P(x) (left) and power spectral density S(f) (right) for the stochastic 
process defined by the stochastic differential equation (|44|) . Dashed green lines are analytical expression ([9]) for the steady-state 
distribution function Po(x) on the left and the slope 1// on the right. Parameters used are s m i n = f, s max = 10 2 , and a = 1 . 




x f 

FIG. 2: (Color online) Probability distribution function P(x) (left) and power spectral density S(f) (right) for the stochastic 
process defined by the stochastic differential equation (|45[) . Dashed green lines are analytical expression ((9J for the steady-state 
distribution function Po(x) on the left and the slope 1// on the right. Parameters used are T m in = 1, x max = 10 2 , and a — 1 . 



The equation was solved using the variable step of integration, Atk — k 2 /x\, with k -C 1 being a small parameter. 
The steady-state probability distribution function Pq(x) and the power spectral density S(f) are presented in Fig. Q] 
We see a good agreement of the numerical results with the analytical expressions. The 1/ / interval in the power 
spectral density in Fig. [T] is approximately between / m j n w 2 x 1CP 1 and / max ~ 2 x 10 2 . The width of this region is 
much narrower than the width of the region 1 <C 2irf <C 10 6 (£ = 10 2 ) predicted in the previous section. 
When rj = —1/2 and v = 3 then equation (JTJ) is 

2 

dx = -2^rdt + -^=dW (45) 

X Z y/X 

We used reflective boundary conditions at x m j n = 1 and x max = 100. The equation was solved with a constant 
step of integration. The steady-state probability distribution function Pq(x) and the power spectral density S(f) are 
presented in Fig. [5] The 1// interval in the power spectral density in Fig. [21 is approximately between / min w 10~ 6 
and / ma x ~ 2 x 10~ 4 . The width of this region is much narrower than the width of the region 10~ 6 -C 27r/ -C 1 
= 10 2 ) predicted in the previous section. 

Numerical solution of the equations confirms the presence of the frequency region for which the power spectral 
density has l//' 3 dependence. The width of this region can be increased by increasing the ratio between minimum 
and maximum values of the stochastic variable x. In addition, the region in the power spectral density with the 
power-law behavior depends on the exponent rj: if 77 = 1 then this width is zero; the width increases with increasing 
the difference \rj — 1|. However, the estimation of the width of the region, obtained in the previous section, is too 
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broad, the width obtained in numerical solutions is narrower. Such a discrepancy can be explained as the result of 
various approximations, made in the derivation. 

VI. DISCUSSION 

In summary, we derived the behavior of the power spectral density from the nonlinear stochastic differential equation. 
In Refs. [29|, |30( only the values of the exponent of the multiplicative noise rj greater than 1 has been used. Here we 
showed, that it is possible to obtain 1/ noise from the same nonlinear SDE for rj < 1, as well. The analysis reveals 
that the power spectrum may be represented as a sum of the Lorentzian spectra with the coefficients proportional to 
the squared first moments of the stochastic variable evaluated with the appropriate eigenfunctions of the corresponding 
Fokker-Planck equation. Nonlinear SDE, corresponding to a particular case of Eq. ((T|) with -q = 0, i.e., with linear 
noise and non-linear drift, was considered in Ref. [9j. It was found that if the damping is decreasing with increase of 
I a; |, then the solution of such a nonlinear SDE has long correlation time. 

As Eq. (l4"Tj) shows, the shape of the power spectrum depends on the behavior of the eigenfunctions and the 
eigenvalues in terms of the function X 2 D(X), where D(X) is the density of eigenvalues. The SDE © considered in 
this article gives the density of eigenvalues D(X) proportional to 1/vA. One obtains l//' 3 behavior of the power 
spectrum when this function XfD(X) is proportional to A - * 3 for a wide range of eigenvalues A, as is the case for SDE 
([3]). Similar condition has been obtained in Ref. [l4|. 

One of the reasons for the appearance of the 1//^ spectrum is the scaling property of the stochastic differential 
equation |T]): changing the stochastic variable from x to x' — ax changes the time-scale of the equation to t' — a 2 ' 1- ' 7 ^ 
, leaving the form of the equation invariant. From this property it follows that it is possible to eliminate the eigenvalue 
A in Eq. (|13| by changing the variable from x to z = X 1 / 2 ^^' ] x. The dependence of the eigenfunction on eigenvalue 
A then enters only via the boundary conditions. Such scaling properties were used estimating the norm of the 
eigenfunction and the first moment X\ of the stochastic variable x evaluated with the A-th eigenfunction. Other 
factor in obtaining the power-law spectrum is wide range of the region of diffusion of the stochastic variable x. 
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